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ABSTRACT 

Dense stellar systems such as globular clusters and dense nuclear clusters are the 
breeding ground of sources of gravitational waves for the advanced detectors LIGO and 
Virgo. The stellar densities reached in these systems lead to the dynamical formation 
of binaries at a rate superior to what one can expect in regions of the galaxy with 
lower densities. Hence, these systems deserve a close study to estimate rates and 
parameter distribution. This is not an easy task, since the evolution of a dense stellar 
cluster involves the integration of N bodies with high resolution in time and space 
and including hard binaries and their encounters and, in the case of gravitational 
waves (GWs), one needs to take into account important relativistic corrections. In this 
work we present the first implementation of the effect of spin in mergers in a direct- 
summation code, NBODY6. We employ non-spinning post-Newtonian corrections to 
the Newtonian accelerations up to 3.5 post-Newtonian (PN) order as well as the spin- 
orbit coupling up to next-to-lowest order and the lowest order spin-spin coupling. 
We integrate spin precession and add a consistent treatment of mergers. We analyse 
the implementation by running a set of two-body experiments and then we run a 
set of 500 simulations of a stellar cluster with a velocity dispersion set to a high 
value to induce relativistic mergers to set a proving ground of the implementation. 
In spite of the large number of mergers in our tests, the application of the algorithm 
is robust. We find in particular the formation of a runaway star whose spin decays 
with the mass it wins, independently of the initial value of the spins of the stars. 
We compare the result with 500 Monte Carlo realisations of the scenario and confirm 
the evolution observed with our direct-summation integrator. More remarkably, the 
subset of compact objects that do not undergo many mergers, and hence represent a 
more realistic system, has a correlation between the final absolute spin and the initial 
choice for the initial distribution, which could provide us with information about the 
evolution of spins in dense clusters once the first detections have started. 

Key words: black hole physics - gravitational waves - methods: TV-body simulations. 



1 INTRODUCTION 

The field of GWs has reached a milestone in the last years 
with the build-up of an international network of GW in- 
terferometers which have achieved their design sensitivity. 
The ground-based detectors LIGO and Virgo are undergo- 
ing major technical upgrades that will increase the volume 



of the observable universe by a factor of a thousand, which 
is referred to as the "advanced" configuration^] 

Dense stellar systems such as globular clusters, galac- 
tic nuclei and, in particular, dense nuclear clusters, are the 
breeding ground of the sources that the advanced detectors 



can expect (see the recent updated review of Benacquista & 
Downing 2011 and also Downing et al.|2011 1. More remark- 



ably, the event rate of stellar-mass black hole binaries, the 
loudest kind of source, will be likely dominated by sources 
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formed dynamically, i.e. via stellar close interactions in these 



stellar systems (Miller & Lauburg|2009 Downing et al. 2010 



Benacquista & Downing 2011 \ 

The data that will be harvested from the advanced 
detectors will allow us to do GW astrophysics. The con- 
struction of templates for matched filtering is crucial in the 
searches for compact binaries. There have been efforts to 
construct these templates by combining post-Newtonian cal- 
culations of the inspiral of the binary with numerical relativ- 
ity simulations of the merger and ringdown. Two appealing 



approaches are the effective-one-body technique ( Buonanno 
& Damour 1999 Buonanno et al.|2009 \ and the phenomeno 



logical hybrid waveform modelling ( |Ajith et al.||2007} |2009| 
Santamaria et al.|2010 |. 

However, the search will be challenging for the simple 
reason that a gravitational wave has not been detected yet. 
Reliable estimates of the event rates for the different kinds 
of binaries and of the expected parameter distribution will 
possibly be crucial for a successful detection. On the other 
hand, once we have the data, we will be able to compare 
the observed rates and parameters with the predictions de- 
rived from different models and thus filter them. This will 
enlighten our understanding of the creation and evolution of 
compact binaries in dense stellar systems. 

The most accurate simulations of dense stellar clusters 
that we can do nowadays are performed with the so-called 
"direct-summation" iV— body algorithms. In particular, the 
family of integrators of Sverre Aarseth has been in develop- 
ment for many decades (von Hoerner 1960 1963| |Aarseth] 
19631. Aarseth's Nbody6 includes both KS regularisation 



(where KS stands for Kustaanheimo-Stiefel) and chain reg- 
ularisation: when particles are tightly bound or their sep- 
aration becomes too small, the system is regularised (see 



Kustaanheimo & Stiefel 1965 Aarseth 2003 \ to avoid too 
small individual time steps and numerical errors. It also 



employs the Ahmad- Cohen neighbour scheme (Ahmad & 
Cohenp973 \ and hierarchical, adaptive time steps. We can 



hence resolve and follow accurately individual orbits in the 
system. In this article we present the first modification of 
a direct-summation code, using Nbody6, that includes all 
non-spinning PN corrections up to 3.5PN order and all spin 
contributions up to 2.5PN order, including spin precession 
equations. 



where v — V1—V2 is the relative velocity vector, m = mi+m-2 
the total mass, r the separation and n = f/r. A and B are 



coefficients that can be found in Blanchet & Iyer ( 2003 ) . The 
spin terms Cn, where N denotes the PN order, are taken 
from ( |Faye et al.|2006[ ) and ( |Tagoshi et al.|2001| >. SO stands 
for spin-orbit and SS for spin-spin coupling. 

These corrections are valid for two isolated bodies and 
shall thus only be applied to the Newtonian acceleration in 
the case of strong, "relativistic" pair-interactions where the 
perturbation by third bodies is sufficiently small. Because 
of this, we deem it reasonable to restrict the implementa- 



tion of PN terms to regularised KS pairs (see |Kustaanheirno 
|fc Stiefel 1965 Aarseth 2003 for details). For this reason 



we also choose the center-of-mass formulation shown in Eq. 
[I] rather than the formulation in the general frame. These 
KS pairs are only formed when the interaction between two 
bodies becomes strong enough so that the pair has to be 
regularised. During the KS regularisation the relative mo- 
tion of the companions is still far from relativistic. Hence, 
only a small, relativistic subset of all regularised KS pairs 
will need post-Newtonian corrections. In order to match the 
order of accuracy of the KS integration in the code, we com- 
pute both the acceleration as shown in Eq. [I] as well as the 
analytical time derivative. 

To save computational costs we switch on the PN cor- 
rections only if one of the following two conditions is fulfilled: 



v > f3c 

v > £c, 
5 



, apN 

and > 7 re i, 

a 



(2) 



where the parameters /3 and 7 are chosen empirically to be 
/3 = 0.02 and 7 = 0.01. Note that this treatment differs from 
Aarseth (20121, who chooses a staggered scheme to switch 



on first PN 2.5, and later PN1 or PN2. We always switch on 
the complete set if equation[2]is fulfilled in order to maintain 
a correct orbit integration under PN influence. 



2.2 Spin Precession 

In addition to the effects on the acceleration, the spin of com- 
pact objects undergoes precession in relativistic two-body 
interactions. This is also taken into account by integrating 
the spin precession equations 



2 THE FORMALISM AND ITS 
IMPLEMENTATION 

2.1 Correction of the accelerations 

We modify the acceleration computation as described in the 



pioneering work of Kupi et al. ( 2006 ) (KAS06) to include rel- 
ativistic corrections, which are based on the post-Newtonian 
(PN) formalism for the interaction between two bodies. We 



note that recently Aarseth (20121 included an approxima- 



tive implementation for relativistic corrections in the new 
version of his code, NB0DY7. The relative acceleration, in the 
center-of-mass form, including all PN corrections used in the 
code can be written in the following way: 



dv 
di 



Gm 



(l+A)n + Bv}+C 1 .5,SO+C2,SS + C 2 .5,SO, (1) 



dS 1 fr 1 t-t 1 t-t 

-jr = -^t/i,so + -gtA.s.ss + -^U2,so 

dE 1 - 1 - 1 - 

-77 = -« Vi.SO + -VV1.5,SS + "tV^.SO 

dt c z c 3 c 4 

S = Si + S 2 

\- I $2 

E = m 1 

ni2 



Si_ 
mi 



(3) 

(4) 
(5) 

(6) 



S and E describe the spin state of the pair. The individual 
terms for Un and Vn, where N denotes the PN order, can 
be found in ( Faye et al.|2006 l and ( Buonanno et al.|2003 l. 



2.3 Relativistic mergers 

Since relativistic binaries lose energy via the dissipative 
2.5PN acceleration term, we need to consistently add a rela- 
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tivistic "merger recipe" in the standard version of the code. 
For the purposes of our study, we must address the following 
points: 

(i) The criterion for two bodies to be transformed into 
one 

(ii) A dynamically correct treatment of the "loss" of one 
star from the simulation 

(iii) Computation of the spin of the star that is formed 
after coalescence from the spins and orbital angular momen- 
tum of the stars that participated in the coalescence 

Post-Newtonian theory can only be applied to the inspi- 
ral of the binary, but not to the actual merger and ringdown. 
We choose for up to 3.5PN order a cut-off distance of 5Rs, 
with Rs = 2G(m\ +m2)/c 2 the combined Schwarzschild ra- 



dius (Yunes & Berti 20081. On the other hand, the newly 



formed compact object must have a mass and a velocity vec- 
tor consistent with the conservation of linear momentum. 
Also, since we are treating spinning compact objects, all 
stars must have an initial spin vector. As we will see ahead, 
in section [3] we use a fitting formula at the last integration 
step before merging the bodies, i.e. at a separation of 5Rs, 
to assign a new spin value to the merged system following 



the prescription of Rezzolla et al. ( 2008 1 



The work we present in this article should be envis- 
aged as a first testing of the algorithm with a "stress test" : 
Our goal is the integration of a large number of relativistic 
mergers in a stellar cluster. We achieve this, as we will see 
later, by setting initially the cluster in a relativistic stage 
with an extremely large central velocity dispersion. In order 
to maximise the number of mergers, we neglect the recoil of 
coalescing pairs, since merging stars with a very large recoil- 
ing velocity could leave the system. However, a priori it is 
straightforward to implement a recipe for the gravitational 
recoil by following a similar fitting formula as in. e.g. the 



work of Pollney et al. (20071; Lousto et al. (20101 



3 TESTING THE IMPLEMENTATION 

In this section we test the implementation itself in a direct- 
summation code. We present tests with a two-body integra- 
tor based on the same routines as Nbody6, but restricted 
to a simple, regularised two-body system. This is exactly 
the part of the modification in the integration that we aim 
at implementing in Nbody6, and hence is a perfect testing 
ground of our algorithm. 

In order to do so, we will compare our simple integra- 
tions with theoretical approaches. In this regard, the for- 



mulae of Peters ( 1964 1 are useful for testing the orbital de- 



cay in the simple non-spinning case. For spinning pairs we 
will check the precession frequencies and conservation of the 
total angular momentum. 



3.1 Non-spinning, merging relativistic binaries 

In this section we compare the results of our approximation 



with the derivation of Peters ( 1964 1 of the evolution of the 



eccentricity and semi-major axis of a binary which is decay- 
ing via the emission of GWs. His derivations are based on 
Keplerian orbits and mimic the 2.5 dissipative term in the 
post-Newtonian expansion. 



Simulation with only 2.5PN 
Peters Formula 



0.45 
0.4 
0.35 
0.3 
0.25 
I 0.2 
0.15 
0.1 
0.05 


1 2 3 4 5 6 

Time (s) 

Figure 1. Comparison of the eccentricity evolution of the two- 
body integration and Peters' approximation. 
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1 + 304 e 



(7) 



In the last equations a is the semi-major axis, e the eccen- 
tricity, t the time, mi and m 2 the mass of the first and 
second star in the binary, G is the gravitational constant 
and c the speed of light. In the case of a circular binary, as 



shown in Peters (19641, one can solve the differential equa- 



tion for a binary with companion masses mi , m 2 and initial 
semi-major axis do: 



1/4 



where 



a(t) = (at - 4/3t) 



64 G 3 mim 2 (mi + m 2 ) 



(8) 



(9) 



This yields a decay time of T c (ao) = ao/(4/3). 

In the general case of eccentric binaries, one can inte- 
grate Eq. |7| numerically and compare the time evolution 
with the results of our simulations. Since Peters' formula 
is only valid for the leading order of gravitational radia- 
tion, we "switch off" the terms 1PN, 2PN, 3PN and 3.5PN 
and only apply the 2.5PN correction. In figure |TJ and |2| 
we show the time evolution of eccentricity and semi-major 
axis for a system with two BHs of masses mi = 10 Mq and 
m 2 = 1 Mq . They agree very well up to the limit of validity 
of the post-Newtonian expansion. 

The 2.5 term only takes into account energy and angular 
momentum loss due to GWs. The 1 and 2 PN terms are 
conservative, they conserve energy and angular momentum, 
and they are the main contribution to periapsis shift. 

In figures |3} and @ we show the time evolution for a 
binary in which we have taken into account the correcting 
terms 1PN, 2PN and 2.5PN. Even though the 1 and 2PN 
terms are conserving energy, the binary coalesces quicker 
than in the Peters approximation, because they change 
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Figure 2. Comparison of the semi-major axis evolution of the 
two-body integration and Peters' approximation. 
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Figure 4. Comparison of the semi-major axis evolution of the 
two-body integration with 1PN, 2PN and 2.5PN terms and Pe- 
ters' approximation. 
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Figure 3. Comparison of the eccentricity evolution of the two- 
body integration with 1PN, 2PN and 2.5PN terms and Peters' 
approximation. 



the orbital velocity and thus the 2.5PN term acts slightly 
stronger. 

The contribution at 3PN and 3.5PN order are small 
compared to the leading order, but these terms cause the 
orbit to diverge when the binary enters the last few 
This is an important effect, since with PN terms up to order 
2.5 one could in principle let the system evolve until an 
overlap of the Schwarzschild radii. When including 3PN and 
3.5PN, on the other hand, this becomes impossible and in 
order to avoid unphysical, divergent behavior one has to 
abort the integration at larger separations. For this reason 
we choose the criterion r = 5Rs where r is the instantaneous 
separation and Rs is the combined Schwarzschild radius. 



Private communication with Seppo Mikkola and Cliff Will 



3.2 Spinning binaries 

3.2.1 Precession of angular momenta 

In post-Newtonian theory, the Newtonian angular momen- 
tum Ljv = x x p, with p = r x mv, is no longer conserved. 
In the case of non-spinning bodies, the direction of Lm is 
conserved and only the modulus Ljv is gradually radiated 
away during inspiral. However, in the case of spinning bod- 
ies this no longer holds (Kidder 1995). Nonetheless, as in 



electromagnetic theory, both the total spin vector 5* and the 
angular momentum vector L precess around the total angu- 
lar momentum vector J — L + S . The angular momentum 
vector we use differs from the usual Newtonian definition: 



L = Ljv + Lipn + Lso + ^2pn. (10) 
With this definition, J = up to 2PN order. The 2.5PN or- 



der, however, introduces radiation loss. Kidder ( 1995 1 esti- 
mated the precession frequency to lowest order, i.e. L = Ljv. 
In the case of a single spinning body with mass m s in a sys- 
tem with total mass m, the precession frequency of both S 
and Ljv is given by 



lc i r i V m a 



(11) 



As an example, let us consider a system of a maximally 
spinning black hole of mass m s = 10 Mq and a non-spinning 
companion of mass 7712 = 1 Mq. We set the system on a 
circular orbit in the x-y plane with radius 10 8 cm with the 
initial spin of m s in x-direction. This gives a total initial 
angular momentum of 



\J\ = y/L z (t = 0) 2 + Si, x (t = 0) 2 = 1.12x10' 



t kg m 2 



(12) 



and thus a precession frequency of ui p — 0.18 Hz. 

From figure ([5| we can see that the approximate value 
for the period of the first precession cycle is (40.4 ± 0.4) s. 
This gives a value of ui p , s i m = 0.15 Hz. The small difference 
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comes from the fact that the calculation assumes the approx- 
imation L — Ljv, and we are already in a very relativistic 
regime. 

Even under the presence of spin-orbitprecession, the 
direction of Jjv should be conserved. Fig. {6} shows the x- 
y-projection of Jjv and Ljv during an inspiral. One can see 
that the direction of Jm is approximately constant but that 
the modulus shrinks due to gravitational radiation. During 
this process, Ljv precesses about this direction. One can also 
see the wobbles in the precession of the orbital plane given 
by Ljv, as described in the appendix of ( Kidder|l995 1. This 
is due to the fact that in reality the corrected L from Eq. 
( 10 1 does the strict precession, which is not true for the 
Newtonian value Ljv, and hence leads Lm to wobble about 
the conserved L. 

The check of J conservation is a powerful way of testing 
the consistency of the approach to estimate the spin and 
angular momentum in the code. 



3.3 Final spin approximation 

In our code we are subject to the limitations of our post- 
Newtonian approach, which is not valid anymore when 
the relative speed becomes larger and larger, i.e. a few 
Schwarzschild radii before the merger. For this, we adopt 



the fitting formula of Rezzolla et al. (2008), derived from 



numerical simulations that address in full general relativity 
the last orbits of the binary, including merger and ringdown. 
We hence implement in the code the following formula for 
the modulus of the final spin (Rezzolla et al. 112008 1 



flfin =73— — s ai + 02 q + 2 a 2 ai \q cos a 

(i + l) 2 (13) 

+ 2(\S 1 \cos P+ \a 2 \q 2 cos j)\l\q + \l\ 2 q 2 ] 1/2 , 

where q = m%lm\ is the mass ratio, Si and S2 the dimen- 
sionless spin vectors and the angles are defined as 



cos ct = 0,1 • 0,2 , 
cos p = &i • I, 

cos 7 = &2 • I. (14) 

Therefore, so as to derive a value for the spin after merger, 
we need the individual spin vectors Si, S2 and the orbital 
angular momentum (OAM) at an arbitrary point in time 
during inspiral. I is a function of the OAM, given by 



|/| = /, o n ( I ^1 1 2 + \a 2 \ 2 q 4, + 2\Si\\a 2 \q 2 cosa) 
(1 + q z Y 

. fss-q + to + 2\ 0,1- 1 2 \ 

I l +q 2 ) (l a l| cos P + \ a 2\q COSJ) 



+ 2V3 + t 2 r) + Urf 



(15) 



where we use the fitting factors , U given in Rezzolla ct al. 
(2008). With Eq. (13} to {15} in hand, one can check whether 
in the regime in which PN is valid, the simulation is consis- 
tent with this formula, in the sense that: 

(i) the total angular momentum must converge to the pre- 
dicted absolute value 

(ii) the predicted final value should be independent of the 
time until coalescence. 

Figure (7} shows the time evolution of both the pre- 
dicted absolute value of the final spin at any given time 
during the inspiral and the actual total angular momentum. 
As one can see, for equal masses this gives a consistent value. 
J is decreasing due to gravitational radiation until it reaches 
the prediction. 



3.4 Energy conservation 

Since Nbody6 is a code to integrate Newtonian systems, it 
regularly checks whether the total energy of the system is 
conserved within some tolerance for numerical errors. In this 
work we have added relativistic terms in the PN approxima- 
tion, so that this is no longer the case: (i) The dissipation, 
mainly by the 2.5PN term, causes a cumulative energy loss 
that has to be tracked and subtracted from the total energy. 
On the other hand, (ii) even the non-dissipative terms cause 
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oscillations in the Newtonian energy, since only the modified 
expression, 

E — i/Ncwt + -E2.5PN.dis + -Eipn + -E2PN + -E3PN + ••• (16) 

is conserved at any given time. We thus calculate and sub- 
tract the corrections up to 3PN order from the total energy 
in order to construct the conserved quantity E. In this way 
we are able to verify energy conservation in the same way 
as it is usually done in purely Newtonian codes. This works 
well if the relativistic corrections are small. However, when 
opn/ci ~ 1 the error induced by PN corrections will domi- 
nate and it becomes impossible to verify the correct integra- 
tion of the system. In order to avoid this, one could decide 
an even larger distance threshold for merging two bodies 
into one or a criterion based on the relative strength of the 
PN corrections. 



4 STELLAR-MASS BINARY MERGERS IN A 
CLUSTER: SOURCES OF GWS FOR 
GROUND-BASED DETECTORS 

It is well-established that most galaxies should harbour a 
massive black hole in their centre, with a mass of some 



1O 6_9 M (see e.g. 


Ferrarese & Ford 


2005 


Ferrarese et al. 


2001 Kormendy & Gebhardt|2001 1. The densities observed 



may even exceed the core density of globular clusters by a 
factor hundred, and hence achieve about 10 7 — lO 8 A/0 pc~ 3 . 
Mass segregation creates a flow of compact objects towards 



2000; Khalisi et al.||2007| |Preto & Amaro-Seoane 2010 


Amaro-Seoane & Preto 


2011 


I, and may build up a cluster 



which could reproduce the effect of an MBH. Indeed, this 
has been used as an alternative to explain phenomena re- 



lated to cluster evolution, like Gl and M15 (Gebhardt et al. 
2002) |Baumgardt et al.|2003a|b||van der Marel et al.|2002[ ). 
Nonetheless, for a globular cluster, compact objects such as 
stellar black holes are very likely expulsed via three body in- 
teractions (Phinney & Sigurdsson|1991 Kulkarni et al. 1993 



0.0!) 




r/r rU (T = 0) 



Figure 8. Mass of the runaway body, M runaw , for each setup, 
averaged over 500 runs. M C \(T = 0) is the total mass of the 
cluster at the time T = and T r i x (T = 0) the initial relaxation 
time of the cluster. The shaded area shows the standard deviation 
for the a = case. 



[20(30l ).[Lee| ( [19951 ) Proved that for a > 300 km/s the merger 
induced by gravity loss in dusters with two components is 
shorter than the required timescale for a third star to inter- 
act with a binary, so that clusters with higher velocity dis- 
persions will not run into that problem. In this section we 
will test the robustness of our code by running simulations 
of dense stellar clusters with a very high velocity dispersion 
to trigger a large number of relativistic coalescences. 



4.1 Initial conditions 

To run a stress test on our implementation, we will consider 
that the clusters are represented by an isotropic Plummer 
sphere containing N — 1, 000 stellar remnants of equal mass 
m. We use N— body units and choose a scaling according to 
KAS06 to trigger a significant amount of relativistic mergers 
to test the code. We set the central velocity dispersion to 
(j ccn w 4, 300 kms -1 , which is equivalent to fixing the ratio 



1 

70' 



(17) 



|Sigurdsson fc Hernquist||1993| |Portegies Zwart fc McMillan| 



In other words, the speed of light "in code units" is c = 70. 
We consider therefore a cluster of compact objects with the 
same mass, spinning with a dimensionless spin parameter 
a and we consider three different initial spin setups for the 
compact objects at the time T = 0: 

• Non-spinning (a = 0) 

• Maximally spinning in the z-direction (a — 1) 

• Random magnitude and orientation 



4.2 Runaway growth 

Because our system consists of very relativistic objects, al- 
most any binary that forms and is regularized will undergo 
a quick merger due to the loss of orbital energy due to the 
dissipative 2.5PN term. Around the time of the core col- 
lapse, i.e. after some ~ 15T r i x (T = 0), with T r i x (T = 0) 
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Figure 9. Spin of the runaway body in each simulation, averaged 
over 500 runs. The shaded area shows the standard deviation for 
the a = case. All initial spin setups lead to a similar evolution, 
except for the very first data point which is slightly higher for the 
maximally spinning initial conditions. 



Figure 10. Cumulative relative energy error in a typical simu- 
lation. In this case we have 71 mergers, which cause the Newto- 
nian energy error to grow significantly. Our alternative method 
to check for energy conservation leads to smaller fluctuations. 



the initial relaxation time of the cluster, a series of mergers 
leads to the formation of one particular star in the system 
that rapidly grows in mass and becomes much more massive 
than the other stars. Therefore we say that the object runs 
away in mass. This is a consequence of the increase in cross 
section for GW capture. The time evolution of the mass of 
this runaway object is shown in figure [8] As we can see, af- 
ter some ~ 15T r i x (T = 0), the runaway object has achieved 
~ 5%M c i(T = 0), a value similar to the case studied in 
KAS06, their figure 1 around 450 time units. 

An important issue that we need to address is the en- 
ergy conservation in the simulations. In figure [10] we show 
both the usual Newtonian energy and the corrected value, 
computed with Eq. 1161. This run includes 71 mergers and 
the Newtonian energy error grows with every single one of 
them. The corrected value for the energy conservation in our 
approach fluctuates significantly less. Nevertheless, we still 
have errors of the order of the total energy in the system. 
We note that in this work we are testing our implementation 
of the relativistic mergers. In a simulation with a realistic 
velocity dispersion, the energy error in coalescences will not 
dominate and hence the check for the conservation of the 
energy in the system would behave much better. 

In order to be able to make a statistical comparison be- 
tween each of the 3 spin setups and the potential impact on 
the evolution of the runaway body, we perform 500 simula- 
tions for each initial spin setup and show the mass averaged 
over each time bin. We can see in figure [9] the evolution of 
the spin for all 3 cases against the accumulated mass of the 
runaway object. Its formation is approximately the same in 
all three different scenarios, and consistent with the results 
of KAS06. Nonetheless, the precise point in time where the 
onset of the runaway process takes places depends sensi- 
tively on the scaling. In any case, the choice for the initial 
distribution of spins is washed out and all three cases show a 
consistent evolution for the runaway body. We additionally 
perform 500 Monte Carlo realisations of the scenario where 
one object merges with non-spinning compact objects com- 



ing in at random angles using the same final spin prediction 
as in the N— body code, so that we can test the statisti- 
cal study. We depict the Monte Carlo spin evolution in fig- 
ure [9] and confirm that this evolution is consistent with our 
N— body analysis within some scattering. 

4.3 Evolution of individual spins 

We now focus on the compact objects that have experienced 
only a few mergers. While the evolution of the spin of the 
runaway object quickly washes out any information regard- 
ing the initial spins, in the case of the other compact objects 
that do not undergo so many mergers, there is a dependence 
on the initial configuration even after core collapse. This is 
particularly interesting, since a trend in the evolution of the 
spin measureable with the advanced detectors would pro- 
vide us with valuable information about the spin evolution 
of compact objects in clusters. 

In figure [Tl] we show the end distribution of spins for 
different initial configurations of the spin distribution for an 
otherwise identical system. The configuration which initially 
had no spins is useful for comparison with the other systems. 
While the x-, y, and z-components individually show no clear 
trend, the absolute value is a a b s = (0.69 ± 0.02). If we move 
on to the second configuration, in which we initially assign 
all compact objects a spin but of random value, the final dis- 
tribution is scattered around the same value, displayed with 
a red line in each of the panels at 0.68. In this case, the final 
value and standard deviation are a a b s = 0.71 ±0.03. Finally, 
if we give all compact objects initially a maximum value 
and set them in a preferred direction, which we arbitrarily 
choose to be the positive z-direction, the final distribution 
has a value of a a b s = 0.76 ± 0.08. 

In figure [12] we can also see this dependence. In the 
plot we display the time evolution of the total spin angular 
momentum in the cluster. In the case of an initially non- 
spinning configuration, the spin builds up from orbital an- 
gular momentum and converges to a generic value in a sim- 
ilar way to what we showed in figure [TT] We are limited in 
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Figure 11. Spin distribution for those objects that have under- 
gone at least one merger during the whole evolution of the cluster 
but not more than four. In the top panels panels we show three 
different initial choices for the spin of the stars. From the left to 
the right we have first a cluster in which initially the stars do not 
have spin, then a random value and in the last column a max- 
imally spinning configuration around the z direction. From the 
top to the bottom panels we display the x-, y-, z- and absolute 
component of the spin ranging between -1 and 1 (a x , a y , a z , a a b s 
shown on the left y-axis of the panels, respectively). The red lines 
depict the values -0.68 and +0.68. 




Figure 12. Dashed green: Total spin angular momentum for a 
cluster in which the remnants are initially maximally spinning in 
the z-direction. Solid blue: Total spin angular momentum for an 
initially non-spinning configuration. 



our analysis to derive the exact value to which the curve 
converges because of an accumulation of numerical errors. 



5 CONCLUSIONS 

In this work we have presented the first implementation of 
the effect of the spin for the treatment of relativistic merg- 
ers in a direct-summation iV— body integrator. For that, we 
modify the calculation of the gravitational forces among par- 



ticles using post-Newtonian up to 3.5 post-Newtonian (PN) 
order and the spin-orbit coupling up to next-to-lowest order 
and the lowest order spin-spin coupling. 

We then check our implementations by running a se- 
ries of tests to compare with results based on analytical 
derivations, for isolated two-body binaries and confirm the 
robustness of our approach. We also present a way to check 
for the correct integration of a system of N particles based 
on tracking the total energy, a usual test with this kind of 
integrators. Our method is valid provided the number of 
relativistic mergers in the system is low. 

The final acid test of the implementation is to compare 
the global dynamical behaviour of a relatively large num- 
ber of stars with the new relativistic behaviour for binaries 
with well-known results based on similar approaches. More 
specifically, we run a similar test to that of KAS06 and ob- 
tain very similar results, which confirms the correct incor- 
poration of the new terms in the code, since the initial spin 
distribution does not significantly change the global evolu- 
tion of the system. This is so, because if two non-spinning, 
equal mass compact objects merge, the merger product will 
be spinning with a w 0.68 ( jDamour fc Nagar||2007| ) in the 
direction of the angular momentum. Since in a Plummer 
sphere there is no preferred direction in the distribution of 
the two-body angular momenta, this leads to a randomisa- 
tion of the non-spinning distribution quickly. In the scenario 
of two maximum spins in the z-direction, i.e. individual spins 
of S = Gm 2 /c with equal masses m, the approximate an- 
gular momentum in the last stable orbit before merger is of 
the same order and thus also rotate the spins and similarly 
wash out the initially preferred direction. 

For the larger subset of compact stars that undergo a 
lower number of coalescences, which is more interesting since 
it is closer to what one could expect to see in a realistic 
cluster, we find that the evolution of the spin for consec- 
utive mergers has a trend that oscillates around the value 



predicted by Damour & Nagar (20071, but with a scatter 



that is a fingerprint of the initial distribution of the isolated 
stars, before they merged with any other in the system. This 
is particularly interesting, since this trend is what will de- 
termine the value of the spin that one can expect to see in 
globular clusters, and should be carefully assessed when de- 
veloping the waveform banks to do the data analysis for the 
first detection. 

Although the systems that we have explored in this 
work cannot be envisaged as representative examples of the 
grounds for which we expect the advanced detectors to ob- 
serve relativistic mergers, the initial study of the behaviour 
of the code is a requirement before we proceed to more real- 
istic systems, and has provided us with initial results which 
could play a crucial role in detection. 

In particular, an immediate goal of our next research 
will be the study of the spin distribution and evolution in a 
dense stellar cluster with a realistic number of stars and in- 
cluding stellar evolution and primordial binaries, such as in 



(Downing et al. 2010 2011 1, but with a more accurate direct 
iV— body integrator. The history and distribution of black 
holes in a dense star cluster is also important for observing 
them in the electromagnetic windows, since it determines 
e.g. number and distribution of X-ray binaries and encoun- 
ters between black holes and other compact objects such as 
neutron stars or white dwarfs. 
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Giersz et al. (20111 clearly show in their (non- 



relativistic) star cluster simulations using the Monte Carlo 
code that quite a few BHs and BH-BH binaries are formed 
and play a role for the dynamics of the central region. The 
presence of BHs may explain the size differences between 
red and blue globular clusters (IDowning 2012 1 and affect 



the number of blue stragglers in a cluster ( Hypki & Giersz 



20131. These papers also discuss that relativistic recoils after 



merger not only important for the gravitational wave signal 
itself, but it is an important ingredient for correct modelling 
of globular clusters. 

The kind of analysis we have presented in this work will 
soon have interesting applications, taking into account that 
the advanced ground-based detectors LIGO and VIRGO will 
have achieved their desing sensitivity as soon as 2016-2017. 
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